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•Vi  :y:o-r\ie.:  Autoregressive  Moving  Average  lime  Series;  Toeplitz  Matrix; 
Modified  Gholesky  Decomposition 

LANGUAGE 

ISO  Fortran 

DESCRIPTION  AND  PURPOSE 

Let  Y(l),  Y (T )  be  a  sample  realization  of  a  mixed  autoregressive 

moving  average  time  series  !Y(t),  t  =  0,  +  1,  ...}  of  order  (p ,q) (denoted 
ARMA ( p ,  q  )  ).  Thus 

p  q 

5  .t(J)Y(t-1)  =  }  B(k)c(t-k)  ,  t  =  0,  +  1,  ... 

j=0  k=0 

for  constants  p,q  ,  u(0)  =  B (0 )  “  1  ,  i(l),  ....  n(p),  B(l),  ...,  8(q)  , 

where  t:  ( • )  is  a  white  noise  series  of  zero  mean  uncorrelated  random  variables 

2  VP  j 

having  variance  o  .  The  zeros  of  the  complex  polynomials  g(z)  =  )  _  pa(j)zJ 

rs  k 

and  h(z)  =  /^_^P(k)z  are  assumed  to  be  outside  the  unit  circle. 

Subroutine  MXPD  calculates  exact,  memory  h,  horizon  t,  minimum  mean 
square  error  linear  predictors  Y(t+h|t)  and  (optionally)  prediction  variances 
j  of  Y(t+h)  given  Y(l),  ....  Y(t)  foi  h  =  h^,  ....  h.,  and  t  *  t^,  ...,  t 
Thu;;  Y(t+h|t)  is  that  linear  combination  of  Y(l),  ...,  Y(t)  closest  to 

i  2 

Y(t+li)  in  the  mean  square  sense  and  ^  =  E{Y(t+h)  -Y(t+h|t)}  is  the 
attained  minimum  mean  square  prediction  variance. 

If  q  =  0,  Y  is  a  pure  autoregressive  process  of  order  p  (AR(p)),  while 
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if  p  =  0  ,  Y  is  a  pure  moving  average  process  of  order  q  (MA(q)).  Subroutines 
ARPD  and  MAPD  calculate  Y(t+h|t)  and  (optionally)  a"  for  the  ARip)  and 
MA(q)  cases  respectively.  Separate  subroutines  are  used  for  these  cases 
due  to  significant  simplifications  in  their  algorithms  from  the  general 
A  RMA ( p , q )  case. 


NUMERICAL  METHOD 

The  method  is  hast'd  on  numerous  special  pr •iportles  of  the  modified 

Cholesky  decomposition  (MCD,  see  Wilkinson  (1967/)  of  the  autocovar l ance 

matrix  of  an  ARMA(p,q)  process.  (See  Newton  and  Pagano  (1930)  for  details), 

Let  [A].j  denote  the  (i,j)th  element  of  t  ,e  .nacr.x  A  and  define 

h  T0EPL{a(0),  ....  a(K-l) }  to  be  the  (K*K)  symmetric  Toeplicz  matrix 

having  =  a ( |  i— j  | )  .  Let  A^,  =  LKDKL^  ,  K  >  i  be  the  MCD  of  the 

symmetric  positive  definite  nested  matrices  A^,  A,  i,_e.  [A^  ^  L 

=  ,  for  i,j  <  K  ,  Lr  is  a  unit  lower-  t:  .an.-,.;.  (K*K)  matrix,  and 

D  is  a  diagonal  matrix.  Then  the  sequences  o;  matrices  I,  .  D 
K-  K '  K.  \ 

and  are  also  nested  and  we  write  (for  example;  ae  (i,j)th  element 


of  any  L^,  for  K  >  max  (i,j)  as  ( I .  ] .  ^  . 


Thus  [D) 


11 


[A] 


11  ’ 


"  'U  ‘  "Al,j 


i-r 

).  f  R  i  j  ^  i  u  i  p  f/  r  l  i  // . 


€-i 


i-l 


[Olii  -  [Aln 


l  lD]ff!l  ■  I'f, 


1  •  K 


We  define  the  following  quantities  for  K  1 

a)  r  -  TOE  PL  (  R  (0 ) ,  -  R..  (K-  1 )  1  -  \...  if 

/. ,  K  /.  /.  /  ,  K  /  .  K  ,  K 


Rz(v)  -  E(  Z ( t ) 7. ( t+v)  ) 


».  +  I.  . 


who  re 

and  ?.(■)  is  an  AR(p) 


-J- 


process  having  coefficients  ot(l),  ...»  u(p)  and  white  noise 
variance  <1  ^  . 


b) 

II 

X  V 

(X(i),  .... 

x(k))t  =  i-;K 

Y  ,  where 

K 

T 

Y  = 
-K 

(Y(l),  .... 

Y  (K) )  . 

c) 

r 

X,K 

-  Elwl>  - 

T 

lx,kdx,k!'x,k  ■ 

Note  that 

d) 

II 

0>> 

(e  ( 1)  ,  .... 

«(k»t  .  k-;K 

4  • 

-1  -T 

=  t  r  t 
X,K  Z,K  Y,K  Z,K 


Then 


Y ( t+h | t )  =  X ( t+h | t )  -  l  m  (  j  )  Y (t+h-j | t) 

j  =  l 


-2  b;1,.  .  ,2 


°  t , h  ^  * LZLX ] t+h , t+h-k  *  °X  ^  t+h-k , t+h-k 


k=0 

".jo  i'V  t+h-k 


where  Y(t+h-j|t)  =  Y(t+h-j)  if  j  >  h,  and 
X ( t+h | t )  = 


,  t+h-k  Jt+h_ktLzl  t+h,£[LX]£,  t+h-kj  | 


q 

X  lLX]t+h,t+h-ke(t+h-k)  ’  h  =  1 . . 

k=h 


,  h  >  q 


Thus 


2 

to  obtain  Y(t+h|t)  and  o^  ^  for  h  =  h^,  ....  h2  and  t  =  t^,  t2 


,  -1 


we  need  LZ,t2+h2  *  4;t2+h2  *  LX,t2+h2  *and  DX,t2+h2  *  Significant  reductions 
in  computing  time  and  storage  requirements  in  obtaining  the  elements  of 
these  matrices  are  afforded  by  noting: 


NOTE __1  Comput  ing  \... 


-1 


Z’  t2+ll2 


The  jth  row  of  b*  is  given  by 

/. ,  K 


'1  = 


( 1  ’  -K-l)  ’ 


....  aj  jd),  1,  Qk_j)  »  2  1  j  i  p 


(0 


»  n(p),  ....  1,  p+1  ±  j  ±  K 


(1) 


(2) 
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where  =  a(^)  if  k  >  p  and 


nj(i) 


«j+,(n  -  ,ij+1  Ci+Ortj^j  d+i-i) 


1 . J 


Tlius  there  are  only  p(p+l)/2  distinct  nonzero, nonone  elements  In  ^ 
K  _>  p  +  1  and  X(l)  =  Y(l)  while 


X(j) 


j;1 

V(j)  +  l  a  Al)'i(l-l) 
1=1  J_1 


j  =  2, 


P 


Y(j) 


+  f  a(l)Y(J-l) 
1=1 


,  j  >  P 


(3) 


NOTE 


2  Computing  h 

•  L  L 


Let  Y  (0 )  =  1,  y(l),  y(2) ,  ...  be  the  coefficients  of  the  MA(«>) 
representation  of  Z ( - )  ,  i.e. 

j'1 

Y  ( j )  =  -  j  a(j-l)yU)  ,  j  >  1  . 

£=max(0, j-p) 


Then 


lLz]p+j,p+j-^  =  yii) 


[LZ]p+)  k  =  -  I  *(r) [L  ] 

Zp+i'k  r-1  Z  P+J-r,£ 


,  0  <  t  <  j  <  K-p 

,  l  -  1 . P-1,  j  >  1 


I...  »  (L*)"1  *  [1 

2 ,  p  Z ,  p 


JZ]  .  J.  ^zJir^r.j-k’ 

%  r=1-V-»-i 


P-» 


lim 


1  im 

K - i  <*> 


l  "'z'k.j 


]-i 

K 


2  =  R7(0 )/a2 


(4) 

(5) 

(6) 

k<j<p 

(7) 

(8) 

(9) 


-  >— 


y  Y (k)  =  R_ (0) la 
k=0  L 


<y(0)  =  a2/  I!  (l-«2(j)) 
/  i=l  J 


By  (5)  and  (6)  the  distinct  elements  of  I,_  ,,  are  contained  in  its 

Z,t„+h. 

T  L  1  T 

first  p  columns,  the  last  of  which  is  (0^  1,  y(l),  ...»  y )) 

Thus  MXPD  finds  the  elements  of  L_  by: 

z,t2+h2 

i)  finding  L  =  (Li1  )_1  ,  ii)  calculating  Y(l) . Y^.-p) 

7. ,  p  Z ,  p  1 

(by  (4))  where  M^-p  is  the  smallest  integer  such  that 


y  Y  (j)  -  R„(O)/02  <  6  ,  (12) 

j=0 

iii)  use  (6)  to  obtain  remaining  elements  of  rows  p+1,  of 

first  p-1  columns  of  L  ,  .  Note  that  (10)  says  that  such  an  M. 

''♦t2+h2  1 

exists  (which  hopefully  is  much  less  than  (8)  and  (9)  say 

that  all  further  elements  of  the  first  p  columns  of  L„  ,,  are  arbitrarily 

Z,t2+h2 

close  to  zero.  Thus  there  are  M^p  elements  to  calculate  and  store  for 

LZ,t2+h2  ‘ 


NOTE  3 


Computing  L 


X,t2+h2  ,  X,t2+h2 


To  compute  and  Dv  .  ,,  we  note  that  (defining  otn(0)  =  1) 

X,t2+h2  X,t2+h2  U 

[  r  ]  =  f  1  a,  .(j-m)  y  o,  .(i-OR-a-B),  i,j>l 

J  I  m=max(l,j-P)  ^  £=max(l,i-p) 


R^di-jl) 


,  i,j  >  p,  |i-j  I  i  q 


,  1  f  j  1  P,  i  >  p,  and 
i  -  j  >  q  or  if 
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R  (v)  =  a  l  8(k)6(k+v)  ,  v  =  0,  ....  q  .  (14) 

X  k=0 

Thus  the  distinct  elements  of  i'^  t  are  R^(0),  ....  R^q)  an<3  t*1e 
elements  in  the  first  p  +  q  rows  and  p  columns.  Further  only  Ry(0),  ...» 
Ry(p+q)  are  required  to  obtain  rx  t  +h  .  These  elements  are  obtained  via 
subroutine  MXCV  by  solving  first  for  v  =  0,  max(p,q) 

f  aUJR^v-j)  =  ?  8 (k) Ry  (v-k)  =  0  ,  v  >  q  (15) 

j=0  -  k=v 

where  Ry£(v)  =  E{ Y (t) c (t+v) }  =  if  v  >  0  ,  where  f>  is  the  Kronecker  delta. 


„  min(v,p) 

Ryc(-v)  =  8(v)o  -  l  a(j)RYE  (-v+j)  ,  v  =  1,  q  .  (16) 

j=l 

Then  Rx(max(p,q)+1) ,  Ry(p+q)  are  obtained  by  (15)  for  v  =  max(p,q)+l, 

....  p+q  . 

To  obtain  L  _  and  D  we  note  that  the  pattern  of  zeros  in 

Aj  ^2  *"”2  a^  1 2  >  n 2 

L„  .  is  the  same  as  that  in  the  lower  triangle  of  r„  and  that  th« 

X,t2+h2  A,t2+i»2 

required  elements  of  T  can  be  calculated  as  needed  (by  (13))  without 

X,  p+q 

storing  them. 

Thus  L„  .  is  calculated  and  stored  in  one  matrix.  To  obtain  the  q 
X,p+q 

nonzero  nonone,  elements  of  the  rows  of  the  rest  of  Lv  ,  we  have: 

a»  1 2  ■f‘”2 


I'x'k.k-j  ■  “<« 


*  j  "  1 »  •  •  *  »  ^ 


Urn  [Dx]k>k 


Thus  rows  p  +  q  +  1,  ...  are  calculated  and  stored  in  a  matrix  having  q 


columns  until  the  elements  of  L  ,  D  have  converged  (at  row  M„  say)  ,  i.e. 

XX  4 

|  [i.xlH  ,  -  «<q-J«)l  '  «.  lloxl„  „  -  »2I  <  J  -  V* . »2-l  <17> 

2 »  2  2 
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iv)  Find  and  1  ,jz  M  +p  1  i j  ’  1  *  -  M]  »  1  1  i  1  P  by  (7),  (6), 

and  (12)  (stored  in  the  integer  HONK  and  (M^*p)  matrix  AK7.)  . 

v)  Find  R  (0),  ....  Rv(q)  by  (14)  via  subroutine  MACV  (stored  in 
X  X 

constant  RXO  and  q-vector  RX) . 

vi)  Find  L„  ,  and  I)v  ,  (stored  in  the  (p+q)*(p+q)  matrix  ALX1 
X,p+q  X,p+q 

and  in  the  first  (p+q)  elements  of  the  f^-vector  DX) . 

vii)  Find  M ,,  and  ( I*x 1 1 j  •  ^Dx^ij  *  1  =  p+q+1*  •••»  M2’ 

j  =  i-q,  ....  i- 1  (stored  in  the  integer  MT WO  and  the  (M^q) 

matrix  ALX2  and  the  rest  of  the  (F^-vector  DX) . 

viii)  Find  X  ,  e  by  (3)  and  (18)  (stored  in  the  t_-vectoisX  and  E)  . 

~ t2  'l2 

ix)  Find  Y(t+h|t)  ,  for  h  =  h^,  ...»  hj  ,  t  =  t^,  ...»  t2  (stored 

in  the  ( 1 2~ t j+1) (hj-hj+D-vectors  YPD  where  Y(t+h|t)  =  YPD((t-t1) 

(h2-h1^-l)  +  (h-h1+l))). 

2 

x)  Find  (optionally)  u  h  which  is  stored  like  YPD  in  a  (t^t^+l) 
(i^-hi+D-vector  PVAR  . 

NOTE  6  Taking  advantage  of  convergence 

To  take  advantage  of  the  convergence  in  1.^,  and  the  user 

specifies  an  absolute  convergence  criterion  6  and  integers  IR0WS1, 

IROWS2  as  the  row  DIMENSION  of  arrays  ALZ  and  ALX2 ,  DX  respectively. 

Thus  IR0WS1  and  IROWS2  must  be  chosen  to  exceed  what  can  reasonably  be 
expected  to  be  M^  and  M2  respectively  or  else  a  nonconvergence  failure 
indicator  will  be  returned  in  ll'AUl.T.  Of  course  if  1ROWSI  or  1  ROWS 2  are 
given  the  value  1 2  t',en  convergence  need  not  be  reached  for  the  algorithm 
to  finis  properly.  Setting  IROWSl  or  IROWS2  smaller  than  tjH^  allows  the 

possibility  of  obtaining  predictors  for  long  time  series  with  a  minimum 
amount  of  required  storage. 


-9- 


Furt her ,  if  the  ^  are  not  to  he  calculated,  the  matrix  ALZ  is  not 
needed  and  l ROWS  1  can  be  set  equal  to  1. 

NOTE  7  Algorithm  for  ARPD 

For  q  =  0  we  have  T„  =  r  ,  p  =  D  =?  L  =  I  D 

Z,K  Y.K  X,K  UZ , K  bX,K  ’  UX,K 

=  K  '  T*1US 


Y  ( t+h  |  t )  «*  -  £  a  ( j  )  Y  ( t+h-j  1 1 ) 

j  =  l 


2  h_1  2 

a.  L  =  l  [L  ] 


t  ,h 


())„] 


k=0 


Z  t+h, t+h-k  1  Z  t+h-k, t+h-k 


2  h_1  ? 

=  o  l  Y  (k)  if  t  >  P 

k=0 


NOTE  8  Algorithm  for  MAPD 


TOEPI 


For  p  -  0  we  have  -  IR  ■»  Lz>|<  -  -  IR  and  1^  - 

(■VO).  •••»  ^(q) .  •••»  0}  ,  (so  that  calculation  of  ALX1  is 


avoided),  Ex  K  =  YK  •  Thus 


Y  ( t  +h 1 1 )  = 


Jh  [Lx]t+h, t+h-k  e(t+h-k) 


t  ,h 


B  (k)e(t+h-k) 

k=h 

0 

h-l 


l  IL 


k=0 


X ' t+h , t+h-k [ DX] t+h-k , t+h-k 


o2  l  B2  (k) 

k=0 


,  t+h  <  Me, 


,  t+h  >  M_ 


,  h  >  q 


,  t  <  M„ 


>  t  >  M 
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STRUCTURE 

SUBROUTINE  MXPD  (NP,  NQ,  ALPHA,  BETA,  S1CSQ,  Y,  I  OPT,  NT],  NT2 ,  NHI ,  NH2 ,  NY PD, 
NPVAR ,  NPPNH2 ,  IROWS1,  IROWS2,  TROWS  I,  DEL,  RYE,  RYEO,  RY ,  RYO ,  IP,  YWK,  RZO, 
RX,  RXO,  ALZI,  ALZ,  ALX1,  ALX2 ,  DX,  MONE,  MTWO,  X,  E,  YPD,  PVAR,  1 FAULT) . 

Formal  parameters 


NP 

Integer 

input : 

order  of  AR  part  of  model 

NQ 

Integer 

input : 

order  of  MA  part  of  model 

ALPHA 

Real  Array  (NP) 

i nput : 

coefficients  of  AR  part  of  model 

BETA 

Real  Array  (NQ) 

i nput : 

coefficients  of  HA  part  of  model  • 

SIGSQ 

Real 

input : 

variance  of  white  noise  in  model 

V 

Real  Array  (NT2) 

input : 

data  vector 

10PT 

Integer 

input : 
1 

0 

option  switch  equal  to: 

if  both  predictors  and  variances 

to  be  calculated 

if  only  predictors  desired 

NT1 

Integer 

i nput : 

tj  (first  memory) 

NT2 

Integer 

input : 

t^  (last  memory) 

NHI 

Integer 

i nput : 

h  (first  horizon) 

NH2 

Integer 

input : 

h  (last  horizon) 

NYPD 

Integer. 

input : 

(NT2-NT1+1) (NH2-NH1+1) 

NPVAR 

Integer 

input : 

same  as  NYPD  if  IOPT  =  1,  >1  if  IOPT=0 

NPPNH2 

Integer 

input : 

NP  +  NH2 

IR0WS1 

Integer 

input : 

row  dimension  of  ALZ  in  calling  pro¬ 
gram  (^NP+2  if  IOPT  =1,  >1  if  IOPT=0) 
see  note  6  above 

I  ROWS 2 

Integer 

input : 

row  dimension  of  ALX2 ,  DX  in  calling 
program  ('NP+NQ-tl)  See  note  6  above. 

' ROWS  3 

Integer 

input : 

row  dimension  of  RY,  IP,  ALZI,  ALXI 
in  calling  program  (>NP+NQ) 

DEL 

Real 

input : 

absolute  convergence  criterion  (see 
(12)  and  (17)). 

RYE 

Rea]  Array  (NQ) 

output : 

Ry, (-1) ,  •••,  Ryc(-q) 

RYEO 

Real 

output : 

^(0) 

RY 

Real  Array  (NP+NQ) 

output : 

Ry(  D . Ry(p+q) 

RYO 

Real 

output : 

Ry(0) 

!  C 

Integer  Array  (IR0WS1) 

work : 

YWK 

Real  Array  (NPPNII2 ) 

work : 

li/.O 

Rea  1 

out  pul  : 

!<.,«>) 

I  X 

Real  Array  (NQ) 

out  put  : 

Rx(l) . R/q) 

RXO 

Real 

out  put : 

Ry(0) 

M.ZI 

Real  Array ( I ROWS 3 , I  ROWS  J) 

out  put  : 

b'.P 

ALZ 

Real  Array  (IR0WS1,  NP) 

output : 

If  IOPT  =  1,  ALZ  contains  [Ly]„, 

)  =  1 . P.  i  =  1 . , 

if  I  OPT  =  0,  ALZ  Is  not  used. 
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AI.X  1 

AI.X2 

Him  1  Army 

Real  Array 

(  1  ROWS  1 ,  1  ROWS'!) 

( 1 KOWS2 , NQ) 

out  pul  : 

output : 

*  X ,  p  tq 

H-x  1  j .  ,  1  =  p+q+i,  . • • *  m2  ’ 

iiX 

Real  Array 

(1  ROWS 2) 

out  put : 

i  ^  1-<1 . 1-1. 

(1>X1,  .  i  =  l.  ....  m2 

MONK 

Integer 

out  pul  : 

M  (see  note  6  above) 

M  rwo 

1 nt oger 

output : 

M  (see  note  6  above) 

X 

Real  Array 

(NT 2)  ■ 

output : 

X(l) . x(t2) 

F 

Real  Array 

(NT2) 

output : 

e(l),  ....  e(t2) 

YPD 

Real  Array 

(NYPD) 

output : 

( (Y ( t+h | t) ,  h  =  h1,  ....  h2)  , 

• 

t  =  tl*  ‘ - • *  t2) 

I’VAR 

Real  Array 

(NPVAR) 

output : 

if  IOPT  =  1,  ((o^  ,  h  =  h  , 

. . . ,  h^) ,  t  t j ,  ....  t2 ) 

if  IOPT  =  0,  PVAR  not  used. 


1  FAULT  Integer 

output:  failure  indicator 

FizhtP*'  Indication  a 

Value  of  1  FAULT 

Meaning 

0. 

no  failure 

1 

NP<1,  NQ'-l,  or  IOPT  not  0  or  1 

2 

NTl-'NP+NQ  or  NT1>NT2 

3 

NH1<1  or  NH1>NH2 

4 

NYPD< (NT2-NT1+1) (NH2-NH1+1)  or  NPVAR<1 
or  IOPT  =  1  and  NPVAR< (NT2-NT1+1) 
(N1I2.-NH1+1) 

5 

IROWS1-NP+2  and  IOPT  =  1  or  IR0WS1<1 
or  NPPNII2- NP+NU2 

6 

I  ROWS2  •" NP4NQ+1 

7 

IROWS3''NP+NO 

a 

SIGSt)  0 

9 

Singular  matrix  in  Subroutine  MXCV 

10 

An  n  (j)  1  (see  (2)  ) 

11 

IOPT  =  1,  1 ROWR1 'NT2+NH2-NP,  and 
convergence  not  reaehed. 

12 

Nonpositive  [  I)  ]  .  .  encountered 

A  11 

13 

I ROWS2^ NT2+NH2-NP-NQ  and  convergence 
not  reached. 

SUBROUTINE  ARPD  (NP, 

ALPHA,  SIGSQ, 

Y,  IOPT, 

Nil,  NT2 ,  NHl ,  NH2 ,  NYPD,  NPPNH2 

YWK ,  GAM,  YPD ,  PVAR, 

I  FAULT) . 

h'nvnal 

param’tars 

NP 

Integer 

input 

order  of  AR  model 

ALPHA 

Real  Array 

(NP) 

i  nput 

coefficients  of  AR  model 

si  esq 

Real 

Input 

variance  of  white  noise 

Y 

Real  Array 

(NT2) 

input 

data  vector 
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10PT 

Integer 

input : 

option  switch  equal  to 

1 

if  both  predictors  and  variances 

to  be  calculated 

0 

If  only  predictors  desired 

NT1 

Integer 

Input : 

tj  (first  memory) 

NT2 

Integer 

input : 

t2  (last  memory) 

NH1 

Integer 

input : 

hj,  (first  horizon) 

Nil  2 

Integer 

input : 

h^  (last  horizon) 

NYPD 

Integer 

input : 

(NT2-NT1+1) (NH2-NH1+1) 

NPPNH2 

Integer 

input : 

NP  +  NH2 

YWK 

Real  Array  (NPPN112) 

workspace : 

GAM 

Real  Array  (NPPNH2) 

output : 

r (1 ) . Y (NP+NH2) 

YPD 

Real  Array  (NYPD) 

output : 

( ( Y ( t  +h | t),  h  =  hj,  ....  h2). 

PVAR 

Real  Array  (NH2) 

output : 

if  I OPT  =  1,  a  ,  h  =  1,  . .  . ,  h 

L  ,  fl 

if  IOPT  =  0,  PVAR  not  used. 

I FAULT 

Integer 

output : 

failure  indicator 

Failure  Indications 

Value  of*  I FAULT 

Meaning 

0 

no  failure 

1 

NP  <  1 

2 

NT1  <  NP  or  NT1  >  NT 2 

3 

Nlll  <  1  or  NH1  >  NH2 

4 

SIGSQ  -  0 

5 

IOPT  not  0  or  1  or  NYPD 

<  (NT2-NT1+1) 

(NH2-NH1+1)  or  NPPNH2 

<  NP+NH2 

SUBROUTINE  MAPD  (NQ, 

BETA,  SIGSQ,  Y 

,  IOPT, 

NT1 ,  NT2 ,  NH1 ,  NH2 , 

NYPD,  NPVAR, 

I  ROWS , 

DEL,  RX,  RXO , 

DX,  ALX,  MTWO, 

E,  YPD, 

PVAR,  I FAULT) 

ormal  parameters 

NQ 

Integer 

input : 

order  of  MA  model 

BETA 

Real  Array 

(NQ) 

input : 

coefficients  of  MA 

model 

SIGSQ 

Real 

input : 

variance  of  white  noise 

Y 

Real  Array 

(NT2; 

input : 

data  vector 

(OPT 

Integer 

input : 

option  switch  equal 

to: 

1 

if  both  predictors 

and  variances 

to  be  calculated 

0 

il  only  predictors 

desired 

NT! 

Integer 

input : 

tj  (first  memory) 

NT  2 

Integer 

input : 

t2  (last  memory) 

Nlll 

Integer 

Input : 

h  j  ( first  hot .  zon) 

NH2 

Integer 

input : 

h2  (last  horizon) 

Nlll 

NH2 
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NY  PI) 

Integer 

input  : 

ril  .'-N  i  l  +  li  (  NI12-NII1+1 ) 

NPVAR 

Integer 

i nput : 

Same  as  NYPD  i£  IOPT  =  1, 

■  1  if  IOPT  -  0. 

I  RoWi> 

I nteger 

input : 

row  dimension  of  AEX  and  I)X  in 
cal  1 f ug  program  (see  notes  6  and 

H  above) 

DEL 

Rea  1 

input  : 

absolute  convergence  criterion 

RX 

Real  Array 

(NQ) 

out  put  : 

R .,  ( 1  ) . Ry(q) 

RXO 

Real 

out  put  : 

R  v  ( 0 ) 

nx 

Rea]  Array 

(1  ROWS) 

out  put  : 

i  =  i . m2 

A  l,X 

Real  Arrav 

( L  ROWS , NQ) 

out  put  : 

I'x’ii’  1  =  l’  •••’  M2 ’ 
j  =  i-q,  . . . ,  i-1 

MTWO 

I nteger 

out  put : 

M  ,  (see  notes  6  and  8  above) 

f; 

Real  Arrav 

(NT2) 

nut  put : 

e  ( l ) . e  ( 1 2 ) 

YPD 

Real  Array 

(NYPD) 

output : 

( ( Y  ( t+h  1 1 ) ,  h  =  h1 . h2), 

t  -  t ^ 1 2 ) 

PVAR 

Real  Array 

(NPVAR) 

output : 

if  IOPT  =  1,  ((o£  h),  h  =  hL . 

b^),  t  =  t^,  t2),  if  IOPT  =  0 

PVAR  not  used 

1  FAULT 

Integer 

output ; 

failure  indicator 

Fiilure  Indications 

Value  of  IFAULT  Meaning 


0 

1 

2 

3 

4 


5 

6 

7 

8 


no  failure 

NQ  •"  1  or  IOPT  not  0  or  1 
NT1  "  NQ+1  or  NT1  ">  NT2 
NH1  '  1  or  NH1  >  NH2 
NY PI)  <  (NT2-NT1+1) (NH2-NH1+1) 
or  NPVAR  <  1  or  IOPT  =  1 
and  NPVAR  <  (NT2-NT1+1) (NH2-NH1+1) 

I  ROWS  <  2 
SIGSQ  <  0 

nonpositivo  [!)  ]..  encountered 

\  IX 

IROWS  <  MT2+NH2  and  convergence  not  reached 


A : ix i Lint'll  a  l.<jorithms 

SUBROUTINE  MACV  (NO,  BETA,  SIGSQ,  RX,  RXO,  I  FAULT) 

mil  pnrame  tern 


NQ 

III  "PA 
SIGSQ 


I nt cger 

Real  Array  (N()) 
Real 


Input:  order  of  MA  model 

input  :  coefficients  of  MA  model 

input:  variance  of  white  noise 
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RX  Real  Array  (NQ) 

RXO  Real 

IFAULT  Integer 


output:  Rx0) . Ry(NQ) 

output:  R^,(0) 

output:  failure  indicator  equal  to: 
0  if  no  failure 
1  if  NQ  <  1 


SUBROUTINE  MXCV  (NP,  NQ,  M,  IROWS,  ALPHA,  BETA,  SIGSQ,  RYE,  RYEO,  WKM, 
IP,  RY,  RYO,  IFAULT) 


Formal 

parameters 

NP 

Integer 

input:  order  of  AR  part  of  model 

NQ 

Integer 

input:  order  of  MA  part  of  model 

M 

Integer 

input:  highest  lag  to  calculate 

(M  >  max  (NP.NQ) ) 

IROWS 

Integer 

input:  row  dimension  of  WKM,  IP 

in  calling  program 

(IROWS  >  max  (NP,NQ) ) 

ALPHA 

Real  Array 

(NP) 

input:  coefficients  of  AR  part  of  model 

BETA 

Real  Array 

(NQ) 

input:  coefficients  of  MA  part  of  model 

SIGSQ 

Real 

input:  variance  of  white  noise 

RYE 

Real  AYray 

(NQ) 

output:  Ryc (-1) . R^f-NQ) 

RYEO 

Real 

output:  RyE(0) 

WKM 

Real  Array 

(IROWS, IROWS)  workspace: 

IP 

Integer  Array  (IROWS) 

workspace : 

RY 

Real  Array 

(M) 

output:  Ryd) ,  ....  Ry(M) 

RYO 

Real 

output:  Ry(0) 

IFAULT 

Integer 

output:  failure  indicator 

Failure 

Indications 

Value  of  IFAULT 


Meaning 


0 

1 

2 

3 

4 


no  failure 

NP  <  1  or  NQ  <  1 

M  '  max  (NP,NQ) 

IROWS  <  max  (NP.NQ) 
singular  matrix  encountered 


The  subroutines  DECOMP  and  SOLV  as  described  by  Moler  (1972)  are  called 
by  subroutine  MXCV. 
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RESTRICTIONS,  TIME,  STOKAOE 

Tf  the  zeros  of  g(z)  are  not  outside  the  unit  circle,  then  one  of 
the  nt .  ( j )  will  be  greater  than  or  equal  to  one  in  magnitude  thus  giving 
IFAIILT  =  10  in  MXPD.  If  the  zeros  of  h(z)  are  not  outside  the  unit  circle 
then  an  element  of  will  become  nonpositive  thus  giving  IFAULT  =  12  in 
MXPD  or  I  FAULT  =  7  in  MAPD. 

The  bulk  of  storage  and  computing  time  in  MXPD  is  devoted  to  the 
'p  matrix  and  the  *  q  matrix  1.^.  The  values  and  increase 
as  the  smallest  zeros  of  g(z)  and  h(z)  approach  the  unit  circle. 


REFERENCES 

Moler,  C.B.  (1972).  Algorithm  423.  Linear  equation  solver.  Comm.  Ass. 

Comp.  MacJi,  274. 

Newton,  H.J.  and  Pagano,  M.  (1980).  The  finite  memory  prediction  of 
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Wilkinson,  J.H.  (1967).  The  solution  of  ill-conditioned  linear  equations", 
in  Mathematical  Methods  for  Digital  Computers  II.  A.  Ralston  and  H.S. 
Wilf,  eds,  65-93. 


SUUROUT  1NE  MXPDI  NP , NO.  ALPHA  .BE  T  A*  S  I GS (3.  V.  I  OPT .  N T 1  ■  NT2  .NHI  ,NH2. 
1 NYPO.NPVAR.NPPNH2.  I ROM  SI  . 1  ROMS 2. I  ROM S  3  .  OEL  . R  YE  .  R  YE 0.  R  V .  R  V 0 .  IP . 
I  YMK .RZO.RX.RXO.ALZ I .ALZ. ALX I.ALX2 . OX . HONE . NT WO . X .£ • VPD. P VAR . 
t  IF  AULT  I 

C 

C  THIS  SUBROUTINE  CALCULATES  PREDICTORS  YPO  AND  (OPTIONALLY) 

C  PREDICTION  VARIANCES  PVAR  FOR  HORIZONS  NHI.....NH2  EACH  FOR 
C  MEMORIES  NTI.....NT2. 

C 

DIMENSION  ALPHA! NP ) .BE  TAINQ) . Y( NT  2) .R YE  I  NO  I  »RY(  IR0MS3  )  . 

1  IP(  I  ROM S3  I  . ALZI I  I  ROMS 3 , I  RUM S3) . ALZ I  (ROM  SI . NP ) , Y MK( NPPNH2) • 

1 ALX I  I IROMS3. IR0MS3 ) . A L X2 I  I ROMS2 .N Q I .0X1  I ROWS 2 ) • X I NT2 ) .RX(NQ) . 

I E I  NT 2  I  . YPDI NYPD) ,P VARINPVAR ) 

DATA  ZERO. ONE. EPS/ 0.0.  1*0.1 .E— I  0/ 

C 

C  CHECK  INPUT  PARAMETERS  S 
C 

I F  A  UL  T  =  l 

IF ( NP.LT. I . OR. NO. LT. 1 . OR.IOPT.LT. 0. OR. lOPT.GT. 1 )  GO  TO  999 
If  A  UL  T  =  2 

IF  I  NT  1 .LE.NPFNQ.OR .NT2.L I.NTl )  GO  TO  999 
IF  AULT  =3 

IF( NHI.LT.I. OR • NH2 .LT • NH 1 )  GO  TO  999 
IF  AULT  =  4 

NCK=(NT?— NT  I  M  )  4  | NM?-NH I  F  I  I 
IF| NYPD.LT.NCK)  GO  TO  999 
IF  I  NPVAR.LT . 1 )  GO  TO  999 

IF! NPVAR.LT.NCK.AND.IOPT.EQ.il  GO  TO  999 

I F  A  UL  T  =  5 

IF | NPPNH2 .LT .NP.NH2 I  GO  TO  999 
IF<  IROMSI.LT. NPF2. ANO. IQPT.EQ. I  )  GO  TO  999 
IF  I  I ROMSI.LT.l  )  GO  TO  999 
IFAULT  =6 

IF  I  IROMS2.LT. NPFNQ4-1)  GO  TO  999 
IF  AULT  =  7 

IF | IRDMS3.LT .NPFNQ )  GO  TO  999 
IFAULT  =  B 

IF | S1GS0.LE . ZERO)  GO  TO  999 
C 

C  FIND  RYO.RYI 1 ) .. .. .R YINPFNO)  Z 

C 

NPPNQ=NPFNO 

CALL  MXCVI NP.NQ.NPPNQ. I ROMS 3. ALPH A. BE T A • S 1 GSQ. H YE • R YEO. 

I ALX I . IP.RY.RYO.  IF!  ) 

I F  A  UL  T  =9 

IFIIFI.F0.4)  GO  TO  999 

r 

C  FIND  ALZI.RZO  (ALZI  INITIALIZED.  ALPHA(J.I)  FORMED  IN  ALXIIJ.I). 
C  J.I.LE.NP.  RZO  FORMED.  ALZI  FORMED  FROM  ELEMENTS  IN  ALXI)  Z 

C 

ifault=io 

DU  20  I=I.NPPNQ 
DD  10  J=I,NPPNQ 
10  ALZ  II  I . J)  =  ZERO 
20  ALZ I  I  I . I ) =ONE 
DO  30  I  -  1  .  NP 
30  ALX  MNP,  I)  .ALPHA!  I  ) 

C 

IF  I NP. F 0. I )  GO  TO  SO 
NP  M I -NP- 1 
DO  40  J=1,NPM1 


J  J  =  NPM I  - J» l 
JJP I* JJ* 1 

PART=ALX1 ( JJPI  * JJP  I  1 
IF(PART.GE.UNE)  GO  TO  999 
OEN  =ONE-PAR  T  *P AR  T 
OO  *0  I =1 . JJ 

JJP1MI=J joi-l 

40  AL  XlIJJ.II-I AL  XI  I JJPI • I >  — P  APT •  AL  X  U  J  JP  I  »  J  J  P I  M I  )  1/ OEN 
C 

50  R/O-SIGSQ 

OO  60  J=1.NP 

60  RZO=RZO/(ONE-ALX1 < J. J  >  AALXI  ( J. Jl ) 

C 

IF(NP.FO.l)  GO  TO  80 
DO  70  J=2.NP 
JM1=J-| 

OO  70 

J J=JMI-I Fl 

70  ALZU  J.  ll*ALXl<  J«1  .  JJ) 

80  CONTINUE 
NPP 1 =NPF| 

OO  90  I=NPPl.NPPNO 
I FS  T  = I  —  NPP 1 

DO  90  J=1 .NP 

j j=  ifstfj 

Jl=NP-JFl 

90  ALZU  I  .JJ»  =  ALPHA{  Jl  I 

C 

C  IFIIUPT.EQ.I),  FIND  MONE  ANO  ALZ  (INVERT  ALZ I .  FIND  ELEMENTS 
C  OF  NP  TM  COLUMN  OF  ALZ  UNTIL  CONVERGENCE. 

C  THEN  FILL  IN  REST  OF  ALZI  S 
C 

IF( (OPT.EO.O I  GO  TO  230 
OO  no  I  =  I,NP 
00  100  J=1,NP 

IOO  ALZI I . JI=ZERO 
110  AL Z (  I  . I ) =QNE 

IF  (  NP  *EQ.  t  )  GO  TO  140 
OO  I  JO  J*2.NP 
JM  I  — J— I 

OO  130  K=1,JM1 
C=ZERO 
JMK= J-K 
JMKP|:jMKtt 

DO  120  IR=JMKPI,J 

120  C=C-ALZ( J. IR  J4ALZI < IR  .  JMK ) 

130  ALZ(J.JMK1=C 
140  CONTINUE 
C 

CK=R20/SIGSQ 
NPP  IsNPH 

AL  Z ( NPP I . NP ) s  — ALPHA! I  I 

SUM SO* ONE ♦ ALPHA ( I  I • ALPHA ( I  I 

MUP=  M I  NO I IROVS1-NP.NT2FNH2-NPI 

00  160  J=2.MUP 

LLOWsMAXOIO  .  J-NPIM 

LUP=  J 

ALD  =  ZEHO 

CC=ONE 

00  150  LL*LLOV . LUP 


L*LL- I 


no  n  n 


i 


JML - J-L 
L I =  NP»L 

IFIL.GT.O)  CC  =  ALZ(L  l,NP» 

150  AL 0= ALD— ALPHA (  J ML  I *CC 

NPP  JiKIP*  J 
AL7(NPPJ.NP)=AL0 
SUMSQ=SUMSQ»ALO*ALD 

IF ( ABSISUNSQ-CKl.LT.OCLJ  GO  TO  130 
160  CON  T 1 NUF 

IF! MUP.LT.NT2FNH2-NP)  GO  TO  170 
MUNE -NT  2 FNH2 
GO  TO  190 
170  CONTINUE 
1FAULT=1 1 
MONE  = IROMS1 
GO  TO  999 
180  MONF=JFNP 
190  CONTINUE 
C 

IF(NP.EQ.l)  GO  TO  220 
NPM 1 -NP- 1 
Ml MNP=MONE-NP 
00  210  J  =  1 • Ml  MNP 
NPP  J— NPf-  J 

DO  210  L=t.NPMI 
C=ZCRO 

00  200  IR=l,NP 
NPPJMB=NPPJ- I R 

200  C=C-ALPHA( IR >*ALZ( NPPJMR.L) 

210  ALZINPPJ.L»=C 
220  CONTINUE 
C 

FIND  RXO.RXI 1 ) . . . . ,RX (NO »  : 

CALL  MACVl NQ . BE  T  A . S I GSQ • RX  »RXO *  IF  1  > 

FIND  ALXI.OXI 1 .OXINPFNOI  (NOTE  THAT  ALPHAII.J)  FOR 
C  I  .LE. J.LE. I.LT.NP  (S  IN  ALZI I  I ♦ 1 » I- J*  1  II  : 

C 

230  CONTINUE 

1)0  250  1  =  1  «  NPPNQ 

DO  2*0  J=l, NPPNQ 
240  ALXll I  .J)=/ERO 
250  ALX 1  I  I . I l=ONE 
IF  A  UL  T  =  1  2 
C 

()X<  1  I  =  RYO 

0(1  340  t  =  2.  NPPNQ 

IMl =1-1 

on  290  Jal.lMl 
C  *  ZERO 

00  260  L  =  I  *  I 
OO  260  M=1.J 
(LMM=I ABSIL— Ml 

IF  (  ILMM.EQ.O  1  C<CMLZI(I  •UTALZII  J.M)*RVO 
IF ( (LMM.GT.O)  C=C» ALZ I ( I ,L> *ALZ II J«M|*RVt ILMMl 
260  CONTINUE 

IF( J.EQ.l )  GO  TO  280 
JM I  *  J— I 

OO  270  L*I.JM1 

C  =C— ALXI ( I ,L )«OX(L )* ALX1 ( J,L I 


270 


280  ALXI ( I . J ) =C/DX< J) 

290  CONTINUE 

C=RX0 

IFll.GT.NP)  GO  TO  320 
C=ZFRO 

DO  310  M=1 , 1 
CI-*LZI( I ,M) 

C2=ZER0 

OU  300  L  =  l.l 
ILMM=I ABStM-L) 

IF( ILMM.EQ.O)  C  2  =  C  2F  ALZ I  (  I.L)*HYO 
IFCILMM.GT.O)  C2=C2*ALZII I. L)*RY( ILMM I 
300  CONTINUE 

310  C  =  C  FC I *C  2 

320  DO  330  L=l.IMI 

330  C=C-OX(L  I  *  ALXI  C  1  «L  I  4ALX  I  I  t  *L  ) 

IF(C.LT.EPS)  GO  TO  999 
340  OXIII^C 
C 

C  FINO  MT WO •  ALX2.  AND  THE  REST  OF  DX  (IN  THIS  SECTION* 

C  I  AND  . . I-l  REPRESENT  THE  INDICES  OF  THE  NONZERO* 

C  NONONE  ELEMENTS  OF  THE  MATRIX  LSUBX.  NOTE  THAT  IF  INTEGERS 
C  M.LE.N.LE.NPFNQ*  THEN  LSUBX(N.M)  IS  STORED  IN  LSUBXIN.M) 

C  WHILE  IF  N.GT.NP+NO.  THEN  LSUBXIN.M)  IS  IN  AL X2I N-NP-NO . 

C  NQ-(N-M)-M).  M  =  N-NQ.  ...  .  N-  1)  : 

C 

IFAULT=I2 

I UP=M I  NO! IROWS2-NP-NQ. NT 24NH2-NP- NO ) 

DO  400  I  1  =  1.  IMP 
I =NPPNQ ♦ I  I 
I MNQ= I — NQ 

ALX2I  I  I  .  1  l=RXI  NOI/OXI  IMNO) 

IF(NQ.EQ.l)  GO  TO  370 
DO  360  J  J -2 • NQ 
J= IMNO* J J— I 
NQ1 =NQ- J J*l 
Jl = J-NPPNQ 
C=BXI NQI > 

JJMl =JJ- 1 

DO  350  LL— t  *  J  JM 1 
L= IMNQFLL-I 
J2=NQ-( J-LM-1 

IF1 J.LE.NPPNO)  C I = AL X  I ( J . L ) 

IF! J.GT. NPPNO)  C  1  =  ALX2< Jl . J2) 

350  C=C-ALX2(  I  I.LLI40XID4CI 

360  ALX2I  I  I . JJ)=C/DX( J ) 

370  C=RXO 

DO  3B0  L= I *NQ 
LL* I-NQ4L-I 

380  C=C-ALX2( II *L ) • ALX2 ( 1 t.L)«DX!LL> 

OXII)=C 

IFI 0X1  I I.LT.EPS)  GO  10  999 

C 

DO  390  J J* 1 » NQ 
N0MJP1=N0-J J4I 

390  IF  I AOSI ALX21  I  I . JJl-BETAl NQMJP1  > I.GT.OEL)  GO  TO  400 

IFI ARS(OX(  I  1 -SIGSQ ) .GT.DEL )  GO  TO  400 
MT*0= I  I 4NPPNQ 
GO  TO  420 
400  CONTINUE 

IFI  I  UP  *L  T*  NT  24  NH2- NPPNO)  GO  TO  410 


n  n 


mtwo=nt2fnh2 

GO  TO  420 
0  CONTINUE 
IF AULT  = 13 
MT¥0= IHOWS2-NPPNQ 
GO  TO  999 
J  CONTINUE 

FIND  XIII . XI NT2) .El  1  >....  .EINT2) 

NPP I -NP ♦ I 
XI  |  )  =  V|  1  | 

IFINP.EO.il  GO  TO  450 
HO  440  J=2 • NP 
C-YI J» 

JMl-J-l 

DO  430  L=1«JM1 
JML  = J-L 

>  C=CFALZII  J.  JML)  *VI  JML  ) 

>  xi  j )=C 

>  OO  4/0  J=NPPI.NT2 
C=V( Jl 

00  460  L-l.NP 
JML=J-L 

)  C=C  ALPHA  I  L  )  *Vl  JML  ) 

1  X(J)=C 

El  l  )  =  XC  I  ) 

OO  490  J=2 • NPPNQ 
C=X I J  > 

JM1 =J-1 

OO  A80  LM.JMI 

>  C=C-ALXl< J.L)*E(L) 

)  FIJI =C 

MI.OW-NPPNO*  1 
MUP=M|NOINT2.MTWOI 
OO  510  J*MLOW.MUP 
JJ= J-NPPNQ 
C=X I J ) 

OO  500  L=l .NO 
LL=J-NQFL-l 

>  C  =C— ALX2 IJJ.LI4E  ILL  I 

>  EIJ)=C 

IF( MUP.EQ.NT2)  GO  TO  540 
MUPP 1 =  MUPF I 
OO  530  J*MUPP1.NT2 
C=Xl J) 

OO  520  L=l .NQ 
JML=J-L 

C=C-BETACLI4EI JML) 

>  E(J)=C 

i  CONTINUE 


ino  vpo  : 

NP0PTsNH2-NHI4l 
00  630  NT  =  N  T |  • NT  2 
NSUF AM= I  NT -NT I  )4NP0PT 
NT  MNP  =  N  T  — NP 

OO  550  l  =  I  • NP 
I  I =NTMNP4I 


550 


T 

-  \ 

YXKC  I  I =Y<  I  I  » 

DO  610  NH  =  1  • NM2 
NPPNH=NP*NH 
NTPNH=NT ♦ NH 
CROWLX=NTPNM-NPPNQ 
XTPH=ZERO 

IFINH.GT « NO )  GO  TO  590 
IFI  NTPNH-Gr.MTurO)  GO  TO  570 
00  56  0  K=NH.NQ 
INDL=NQ-K>  l 
INOE=NTPNH-K 

X  TPH=  XTPH*ALX2( I ROWL  X • INDL I *E (  INDE ) 

GO  TO  590 

DO  580  K=NH.NQ 
I NOE=NTPNH— K 

X  TPH- X TPH  +  BE TA(K )*E( INDE) 

C*XTPH 

DO  600  J=l.NP 
I  NO Y- NPPNH—  J 
C=C-ALPH*( J ) *YWK< INOY  | 

Y»K|NPPNH) =C 
DO  620  l=NHl.NH2 
N 1 0=NP ♦  I 

N 1  l - I -NH  t  ♦ t 
N1 2=NS0F ARANI 1 
VPD<N12) =Y  VK(NIO) 

630  CONTINUE 
C 

C  IF  IOPT.EQ.1.  RIND  PVAR  i 
C 

IF ( IOPT.EQ. 0 )  GO  TO  690 
Ml MNP=MONE-NP 
DO  680  NT -NT  I  • NT2 
NSOFAR=f NT-NT1 )*NPDPT 
OO  680  NM=NHI.NH2 
NINDX=NSQFAR*NH-NHi ♦  ! 

NTPNH=NT  ANH 
C=ZERO 

DO  670  KP1=I.NH 
K  =KP 1—1 
KROW=NTPNH-K 
C 1=SIGSQ 

IFIKROW.LT.MTWOI  CI=OXCXROW) 

C2=ONE 

IFIK.EO.OI  GO  TO  670 
nppk=np»k 

C23ALZ<NPPK,NP» 

OO  660  I R* 1 • K 
K  MR  =  K - I R 
C3=ZEPO 

IF ( KMR.GT.M1MNPI  GO  TO  640 
!NDL=NPfKMR 
C  3=  AL  Z < INOL.NPI 
640  C4=ZEPO 

IF(IR.GT.NO)  GO  TO  660 
N10=Nf PNH-KMR 
IFIN10.GT.MTM0)  GO  TO  650 
NI I =NI 0-NPPNQ 
N1 2=NO-IN»l 
C  4  =  AL  X2 ( N 1 1.N12) 

GO  TO  660 


560 

5  70 

580 

590 

600 

610 

620 


6S0  C 4  - OE  T  A  ( I H ) 

a  f>o  c?-C2*rj*c4 

A  70  C  =  C  FC 1 FC  2*C  2 

680  PV  AH  <  N I  NOX  I  =  C 
C 

690  CONTINUE 

c 

ir  aul  t=o 

999  HE TURN 
END 

SUBROUTINE  M  APD  <  NO  , BETA. SIGSQ.Y. IOPT. NTI.NT2.NH1.NH2, 

I  NYPD.NPVAR,  IHOMS.0EL.RX.RXO  rOX  •  AL  X  .  MT  Ml)  .E  •  YPD.  P  VAR  •  I F  AULT  ) 

C 


C  THIS  SUBROUTINE  CALCULATES  PREDICTORS  YPD  AND 
C  (OPTIONALLY!  PREDICTION  VARIANCES  PVAR  FOR  HOR12INS 
C  NHI.....NH2  EACH  FOR  MEMORIES  NT1.....NT2. 


C 

DIMENSION  BETA(NQ) . Y(NT2 ) .RX(NQ) . DXI  I ROMS) . AL XI  I RO MS* NQ I • 
!F( NT2 I . YPO< NY PD) .PVAR ( NP VAR ) 

OAT A  ZERO. ONE. EPS/0.0. 1*0.1 *E- l 0/ 

C 

IFAULT=I 

IF< NQ.LT . 1 .OR. IOPT .Lf. O.OR. IOPT.GT. 1 )  GO  TO  199 
IF  AULT =  2 

IF  (  NT  I  .LT.NQM  .OR.NTI  .GT.NT2)  GO  TO  199 
IF AULT=3 

IF(NH|*LT«I« OR • NH I «GT • NH2 )  GO  TO  199 
IFAULT=4 

NC  K  =  (  NT  2—  N  T  1  ♦'I  )*(NH2-NHI*1  ) 

IF(NVPO.LT.NCK)  GO  TO  199 
IFINPVAR.LT  *  I )  GO  TO  199 

IF < NPVAR.LT.NCK. AND. IOPT .EQ. I )  GO  TO  199 
IF A  ULT  =5 

IF ( IRUWS.LT .2)  GO  TO  199 
IF  AULT=6 

IF< SIGSQ.LE .ZERO)  GO  TO  199 
IF  AUL  T=  7 
C 

C  FIND  MTMO  ANO  ROMS  OF  AL X ,  DX  I 
C 

CALL  MACVINO.BETA. S I GS Q. RX • RXO • IF  I  I 

C 

OxC 1 )=RXO 
DO  10  1=1. NO 
10  AL  X (  I . I )=ONE 
E  (  1  1  =  Y  (  I  ) 

IUP*MIN0I I  ROMS .NT 2 FNH2 I 
DO  100  1=2. I UP 

| | sMAXOC I— NO” 1 •  0  )  F  1 
|MI=I-I 

NEL  TS= I M l -  I l  F  l 

DO  30  J=II.IM1 
J I  NO  = J-  I  I  Fl 
1MJ=I - J 
JMI=J-l 
C=RX<  IMJ  ) 

IF  I J.EO. It )  GO  TO  30 
Jl =MAX0I J-NO-I .  0  )  F  1 
DU  20  L= I  I  <  JM1 
LL=L-I 1 F | 

J J=L  - J I ♦  I 


n  n 


20  C=C-ALX (  I  ,LL  I  *Ox <L >*ALX< J « J J » 

30  AL  X  1  1 • J  I  ND l—C /D X 1 J  1 

C  =  RXO 

on  40  J= 1  1  »  I  Ml 
J  I  N0=  J - I  I FI 

*0  C=C-OX<  J **ALX(  I  ,JIND»**LX( (  .  J( NO) 

IF  1 C.LE .EPS I  GO  TO  199 
OX< I ) =C 


IF  ( 

i.GT. 

NT  2  ) 

GO 

TO 

60 

C=Y 

(  I  1 

DO  50 

J-  1 

.NELTS 

I  1  =  1- 

NFL  T 

SF  J 

-1 

50 

C=C- ALXl  I 

.  J) 

*61 

I  I  ) 

E  (  I 

)  =C 

I  F  ( 

I.LE. 

NQ) 

GO 

TO 

1  00 

60 

OO  70 

J=1 

•  NO 

• 

J  J=NO 

-  JF  1 

70 

IF  1  ABSl  AL 

X(  I 

.  J  ) 

-BETA! JJ | J.GE.OEL 1 

IF!  ABS10X1 I  J-SIGSQ  l.GE.DEL)  GO  TO  100 
Mr WO: I 

IF1I.GE.NT2I  GO  TO  110 
IP  1  =  I  ♦  I 

DO  90  J3 I P l • NT  2 

c=yi ji 

DO  80  K*1.NQ 
JMK=J-K 

80  C=C-BETA|K)*E1 JMK) 

90  E1J1-C 

GO  TO  120 
100  CONTINUE 
I F  A  UL  T  ■=  8 

IF(  IUP.LT.NT2FNH2 1  GO  TO  199 
110  MT  WO— NT  2  FNH2 
120  IFAULT  =  0 

CALCULATE  PREDICTORS  : 

NPDPT-NH2-NH1 F | 

OO  140  NT=NTl  .  NT2 
NSOF  ARs 1 NT -N  T 1 )*NPOPT 
DO  140  NH-NM1.NH2 
N I ND— NSOF  ARFNH—  NH1 F 1 

ypo<nino)=zero 

IF1NH.GT.NQI  GO  TO  140 
NTPNrt-NT FNM 
C=7ERO 

DO  130  K=NH,NQ 
I NDE=NTPNH-K 
I NI)L  =  NQ— KF 1 

C  I  =8E  T  A 1 K  1 

IF( NT  PNH .LE.MTWOI  Cl - ALX  <  NT  PNH* 1 NOL ) 
C-CfC 1 4E 1  INDEI 
TPD INI  NO ) =C 

IF  lUPTxl.  CALCULATE  VARIANCES  1 

IF! IOPT .EQ. 0)  GO  TO  199 
DO  180  NT=NT1,NT2 
NSOF AR=(NT-NT1 ) FNPDPT 
DO  180  NM=NH1.NH2 
NHMI «NM- I 


130 

140 

C 

C 

c 


too 


n  n 


nind-nsofarfnh-nhi  *■  i 
N  T PNM -N  f F NH 

c  =  s  ir.so 

IF  (  NT  »LT  .  M  T  MO  I  C^UX(MPNH| 

IF  (  NH.El).  1  >  GO  TO  t  80 
NHUP-MINOINHMl  .NO  I 
IFlNT.GE.MTMOl  GO  TO  160 
DO  150  K-I.NHUP 
I  NOD*  NTP  NH-K 
I NOL  =  NO—  Kt 1 

l 50  C  =  CFOX< 1NOO I FALX (N1PNH, INOL ) * AL X I NTPNH . I NOL » 

GO  TO  100 

160  0(1  1/0  K=1.NHUP 

1/0  C  =  C  F5 IGSQ*BE  T  A  ( K )• BE  TA(K) 

lflO  PVARININO»=C 
C 

190  PFTUHN 

FND 

SUB ROUT INE  ARPOI NP. ALPHA .SI GSO.Y . I0PT.NTl.NT2.NHl *NH2. 

I NYPO.NPPNH2 . VWK.GAM. YPD .PVAR.IF AUL T  > 

C 

C  THIS  SUBROUTINE  CALCULATES  PREDICTORS  YPO  AND  I  OPT (ON ALLY  I 
C  PREDICTION  VARIANCES  PVAR  FOR  HORIZONS  NH1.....NH2  EACH 
C  FOR  MEMORIES  NTI.....NT2 
C 

DIMENSION  Y(NT2).ALPHA(NP|  . YMK ( NP  PNH2 1 • GAM 1 NPPNH2 )  . 

IYPOI NYPOl .PVAHINH? I 
OATA  ZERO .ONE/ 0 . 0 « 1*0/ 

C 

I F  A  UL  T  =  I 

IFINP.LT.t)  GO  TO  100 
IF AULT  =2 

IF ( NT1  .LE.NP. JR. NT  I .GT.NT2)  GO  TO  100 
IF AULT=J 

IF  1 NH 1  .LT . I • OR .NHl  . GT . NH2 1  GO  TO  100 
(FAULT-* 

IFI SIGSQ.LE.ZERU1  GO  TO  100 
IF  AULT  =5 

IF( I0PT.LT.0. OB. IOPT.GT.il  GO  TO  100 
I F  AOL  T  =  6 

NCK=<  NT2-NT 1  FI  I *( NH2-NH1  F 1  I 

IFI NYPO.LT. NCK.OR.NPPNH2.LT. NPFNH2I  GO  TO  100 
I F  AUL  T  =  0 
C 

FIND  PREDICTIONS  : 

NPDPT  =NH2-NH I F I 
DO  50  NT  =NT 1 . NT  2 
NSOF AR=(NT-NT I 1FNPOPT 
NT  MNP=NT-NP 

OO  10  1=1. NP 
I  I  =  NT-NP ♦ I 
10  YWK(  I  1  =  Y  I  1  I  1 

OO  30  NH=1.NH2 
NPPNM=NP ♦ NH 
NTPNH=NT  fnh 
C*ZERO 

OO  20  I-*l  ,NP 
(  UNPPNH-I 

20  C  =C  — ALPHA (11AYMK1II1 

30  YMK  1 NPPNH I =C 


OQ  *0  NH  =  NM 1 • NM2 
INONH  =  NSOFAR«-NH-NHI M 
I  NDWK-=NP  *■  NH 

40  YPO! INONHI -YWK( INDMK) 

50  CONTINUE 

IF  IOPT  =  l •  FINO  VARIANCES  : 

IF  I  I OPT • EQ. 0 )  GO  TO  100 
GAM (  1  ) =  —  AL  P  H  A I  1  > 
tF<NH2.Ea.ll  GO  TO  80 
OO  TO  NH=2.NH2 
LLO*=MAXO!O.NN-NP1  I 
C=2ERO 

00  60  LL=LL04.NH 
L— LL— 1 
Cl -ONE 

IF(L.GT.O)  C l=G AM(L I 
NHML=NH-L 

60  C=C— ALPHA! NHML I *C 1 

TO  GAM | NH ) =C 
80  PVARl I l-SIGSO 

IF  <  NH2.EQ. I  I  GO  TO  100 
00  90  1=2. NH2 

IM  |  —  I  —  1 

90  PVARl I  1 — PVAR 1  I M 1 |fSIGSO*GAN< I Ml  I* GAM!  I Ml  I 
100  RETURN 
END 

SUBROUTINE  M  At  V  !  NQ  .  BET  A.  SI  GSQ.RY  •  RVO  •  IF  AULT  ) 

THIS  SUBROUTINE  CALCULATES  MAINQ)  AUTOCOVARIANCES  OF  LAGS 
O.....NQ  INQ.GT.01 

DIMENSION  BE TA1 NQ) .RY( NQ) 

DATA  ONE/ l . 0/ 

IFAULT=l 

IFINQ.LT. 1)  GO  TO  40 
IFAUL  T  =  0 

C=ONE 

OO  10  1=1. NQ 
10  C=C*BET A! I 14BETA1  (  ) 

RYO=C*SIGSQ 

DO  30  I  V  =  1 «  NQ 
C  =  BET  A!  IV» 

IF! IV.EQ.NQ)  GO  TO  30 
NQM I V=NU- I  V 

OO  20  J  =  I  • NQM | V 
JPt V= J* l V 

20  C=C*BETA! JI4BETA! JPIV) 

30  RT« IV)=C*SIGSQ 
40  RETURN 
FNO 

SUBROUTINE  MXCVINP .NQ. M, IR 04 S . ALP HA . BE T A , S I GSQ . RYE .R VEO . 
I MKM . IP.RY.RYO. I FAULT) 


THIS  SUBROUTINE  CALCULATES  ARMAINP.NQ)  AUTOCOVARIANCES  FOR 
LAGS  0 . • . • . M  ! M.GE.MAX! NP.NQ) •  NP.NO.GT.O)  i 


DIMENSION  ALPHA! NP ) . BET A !NQ  )  .R YE!  NO I .NKM! IRONS. IRONS)  • 
1 IP< IRONS) ,RV(M) 

OATA  ONF. ZERO/ 1 .0. 0.0/ 

C 

I F  A  UL  T  =  t 

IF< NU.tr. I .OR.NP.LT. 1 )  GO  TO  110 
1FAULT=2 

MAXPQ=MAXO!  NP. NQ) 

MMzMAXPQAl 

IF( M.LE.MAXPQ)  GO  TO  110 
IFAULT=3 

IF( IRONS.LT. MM)  GO  TO  110 
IF  A  UL  T  =4 

C 

C  FIND  RYEO.RYE! I )..... RYE INO)  5 

C 

HYF0=SI GSQ 
OO  30  1V*I.NQ 
C=S IGSQABETA! IV) 

NUP^M INOI IV. NP ) 

DO  20  J= 1 .NUP 
I VMJ-I V- J 

IF! IVMJ.EQ.O)  GO  TO  10 
C=C-ALPHA! J 1 *RY  E I IVMJ  ) 

GO  TO  20 

10  C=C— ALPHA! J)*RYE0 

20  CONTINUE 

30  RYE!IV)=C 

C 

C  USE  OECOMP.  SOLV  TO  OBTAIN  R YO.RY I  1 )..... RY! MAXI NP . NO ) ) 

C 

00  40  I V= 1 . MM 
RY 1  IV  )  =  ZERO 

DO  AO  J - 1 . MM 
AO  NKM( IV. J)=ZERO 

C 

NPP | =  NP  A 1 
NOPl=NO* 1 
DO  60  1 VPl  =  l  *  NQP 1 

I  V  =  I  VP l- 1 
C=RYFO 

irilV.GT.O)  C=CABETA!I V) 

IF! IV.EQ.NQ)  GO  TO  60 
DO  SO  K=IVPI.NQ 
KMI V  =  K-I V 

SO  C -  C  AI1E T  A ! K ) A RYE ! KM I  V I 

60  RY! I VPl )-C 
C 

DU  TO  I  VP 1  =  1  .MM 
I  Vs  IVPl-1 

NKM!  I VPl . I VPl  ) =  NKM!  I  VP  1.  1  VPl |  AONE 
DO  TO  J=I.NP 
I  l-IABSC  IV-JlAl 

TO  NKM! I VPl • t 1  ) =NKMl 1 VPl. I  1 ) A ALPHA I J ) 

C 

CALL  OECOMP! MM, IRONS. NKM.1P) 

IF! IP(MM).EQ.O)  GO  TO  110 
IFAULT-0 

CALL  SOLV ! MM  » I  RONS • NKM  «R V* IP) 

RYO=RV<  I  I 

DO  80  I V*l  «  MAXPQ 


IVP1=IVM 

80  RVI  IV)=Rr(  I VP1  ) 

c 

C  USE  DIFFERENCE  EQUATION  TO  GET  THE  REST  OF  THE  RV  : 
C 

IF ( M.EQ.MAXPQ)  GO  TO  110 

DO  100  1 V=MM . M 

C=2ERO 

00  90  J>1.NP 
I VM J=1V- J 

90  C=C— ALPHA ( J)*RY( IVMJ) 

100  RV (  IV  »=C 
C 


110  RETURN 
END 


END 

DATE 
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